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modified critical scaling in the Kibble-Zurek mechanism of defect formation 


R. G. McDonald and A. S. Bradley 

Department of Physics, QSO — Centre for Quantum Science, 
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As a test of the Kibble-Zurek mechanism (KZM) of defect formation, we simulate the Bose-Einstein con¬ 
densation transition in a toroidally confined Bose gas using the stochastic projected Gross-Pitaevskii equation 
(SPGPE), with and without the energy-damping reservoir interaction. Energy-damping alters the scaling of 
the winding number distribution with the quench time - a departure from the universal KZM theory that relies 
on equilibrium critical exponents. Numerical values are obtained for the correlation-length critical exponent 
V and the dynamical critical exponent z for each variant of reservoir interaction theory. The energy-damping 
reservoir interactions cause significant modification of the dynamical critical exponent of the phase transition, 
whilst preserving the essential KZM critical scaling behavior. Comparison of numerical and analytical two-point 
correlation functions further illustrates the effect of energy damping on the correlation length during freeze out. 


I. INTRODUCTION 

The Kibble-Zurek mechanism (KZM) describes defect for¬ 
mation in the symmetry-breaking dynamics of a system un¬ 
dergoing a second-order phase transition [1, 2]. The theory 
exploits critical slowing down near the transition, whereby the 
system relaxation time diverges and the system becomes es¬ 
sentially frozen, allowing the description of critical dynamics 
in terms of equilibrium critical exponents [3]; thus the quench 
phenomena are set by the universality class of the system. The 
critical exponents determine the density of defects introduced 
by symmetry breaking during the quench, and thus are central 
to the testable predictions of KZM. Bose-Einstein condensa¬ 
tion (BEC) is a t/(l) symmetry breaking transition whereby 
the phase of the order parameter acquires independent values 
over finite-sized domains, the size of which depends on the 
speed of the quench. Theoretical treatments of KZM in Bose- 
Einstein condensation have used a description of the Bose gas 
that can be reduced to a form of stochastic Ginzburg-Landau 
(G-L) theory [4-11], involving a Gross-Pitaevskii equation 
coupled to a grand canonical reservoir providing a source of 
particles and thermal noise. Eundamentally, the model is a 
mean field theory driven by the simplest possible reservoir 
coupling, and analyses of defect formation confirm the equi¬ 
librium scaling hypothesis of KZM. A basic question then 
arises: what is the role of a specific system’s reservoir inter¬ 
actions in determining the critical dynamics during a quench? 

A rigorous and tractable reservoir interaction theory has 
been developed for the dilute Bose gas from first principles 
in the form of the stochastic projected Gross-Pitaevskii equa¬ 
tion (SPGPE) [12-16]. The theory is a synthesis of quantum 
kinetic theory [17] and the projected Gross-Pitaevskii equa¬ 
tion [18], and provides a tractable approach for numerical sim¬ 
ulations of critical dynamics that includes all significant reser¬ 
voir interaction processes. The SPGPE describes the evolu¬ 
tion of a high-temperature partially condensed system within 
a classical field approximation, is valid on either side of the 
critical point, and in 3D has been used to quantitatively model 
the phase transition [8], and high-temperature dynamics [19] 
observed in experiments. The complete SPGPE includes a 
number-damping reservoir interaction (G-L type) described in 


previous works [8, 19], and an additional interaction involv¬ 
ing exchange of energy with the reservoir [20]; the latter is a 
number-conserving interaction that can have a significant in¬ 
fluence on system evolution far from equilibrium [15, 21, 22], 
such as occurs in the region of critical slowing down near the 
transition. 

In this work we investigate the effect of the energy damp¬ 
ing reservoir interaction on the outcome of quenches across 
the Bose-Einstein condensation transition. We consider a 
toroidally trapped Bose gas consisting of a quasi-ID super¬ 
fluid fraction immersed in a 3D thermal cloud. This system 
may be modelled using the effective ID-SPGPE description 
for elongated systems [22]. A finite persistent current circu¬ 
lating around the ring provides a clear signature of symmetry 
breaking [23]. We perform quenches of the reservoir chemical 
potential for a range of quench times, with and without the en¬ 
ergy damping terms, to produce a statistical distribution of the 
final winding number of the persistent current. The Kibble- 
Zurek mechanism predicts a power law relation between the 
quench time and the standard deviation of the final winding 
number. We numerically determine the power-law exponent 
for these relations using the two reservoir interaction theories 
and compare with the predictions of mean field theory. In¬ 
clusion of energy damping causes the power-law exponents to 
depart from the mean field values. In particular, the dynamical 
critical exponent extracted from the complete SPGPE simula¬ 
tions differs from that predicted by mean held theory, does 
not conform to a known universality class, and suggests non- 
universal modifications to the critical dynamics. The effects of 
energy damping are further exemplified by comparing the dy¬ 
namics of condensate number and two-point correlations for 
the two theories. 


II. THEORY 

A. The Kibble-Zurek mechanism 

We first review aspects of KZM relevant for this work. We 
consider a system driven across a second order phase transi¬ 
tion by a quench of the chemical potential from -po to po- We 
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However, as the system approaches the critical point, j-L - 0, 
the relaxation time diverges and there is thus an instant dur¬ 
ing the quench where the relaxation time is equal to the time 
remaining reach the critical point. This time t - -i is called 
the freeze-out time, after which the system cannot keep up 
with the external parameter, and thus remains “frozen” (im¬ 
pulse regime) until time -\-t after the transition when it can 
again respond to the environment, as shown in Fig. 1 (a). The 
freeze-out time satisfies the equation 

t = T(-i), (4) 

which may be solved [23] to give 

f = (To'rg)'"^" (5) 


where 



FIG. 1. (colour online) (a) Adiabatic-impulse model of dynamical 
symmetry breaking in KZM. The boundary between regions t is de¬ 
fined as the instant when the time remaining until the critical point 
is equal to the system relaxation time, i.e. when t = T(-f). (b) 
Schematic of defect formation during Bose-Einstein condensation 
in a ring trap. The scale of regions acquiring a t/(l) symmetry- 
broken phase is set hy the coherence length at the freeze-out time, 
I = Phase differences between adjacent domains, Af?,, are re¬ 

solved through the formation of defects, (c) Reservoir interactions 
in the stochastic projected Gross-Pitaevskii theory of the dilute Bose 
gas. The numher-damping (y) interaction drives condensate growth. 
The number-conserving energy damping process (e) has a significant 
role far from equilibrium. 


zv 

1 -I- zv 


( 6 ) 


The adiabatic-impulse approximation assumes that the sys¬ 
tem ceases to follow the equilibrium solution adiabatically 
precisely at —t, when the dynamics are frozen until the sys¬ 
tem returns to adiabatic following at t. Using the freeze-out 
time we can obtain the freeze-out chemical potential and cor¬ 
relation length 





(7) 


and 


1-^0 




( 8 ) 


The density of topological defects can be estimated as the size 
of the symmetry broken domains at the freeze-out time [see 
Fig. 1 (b)], and is then given by 


define the reduced parameter e(t) 



(9) 



jUO Tq 


( 1 ) 


where 2tq is the quench duration and t 6 \—tq,tq\. The 
equilibrium correlation length and dynamical relaxation time 
are related to the reduced parameter by 


^(0 = 


^0 


and 


T(f) = 


Tq 

k(f)P’ 


( 2 ) 

(3) 


where D is the dimensionality of the system and d is the di¬ 
mensionality of the defects. This is typically an over-estimate 
of what is observed, and thus f is commonly replaced by sf, 
where i ~ (9(1 - 10) depends upon the model [4—6, 10]. From 
(5) it is apparent that tq may be chosen small enough such that 
the freeze-out time is larger than the total ramp time (f > tq). 
This implies a fast quench limit, below which the system will 
effectively experience a jump in the chemical potential, rather 
than a ramp. To be in the regime where KZM scaling should 
apply, the quench time must satisfy tq > tq. 

In a ID Bose gas, the defects that form as a result of merg¬ 
ing domains are grey solitons [9, 24]. According to (9), the 
density of solitons after a quench should scale as 


respectively, where v and z are critical exponents defining the 
universality class of the phase transition, and fo and tq are 
constants that depend on the microscopic details of the sys¬ 
tem. Initially the system follows the quench adiabatically. 


IIsOCTq^^. ( 10 ) 

The problem arises that grey solitons are unstable [9, 25- 
27], and thus in the presence of any dissipation, the solitons 
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will vanish following the quench. This inhibits verification 
of KZM via soliton counting as it is not clear how long af¬ 
ter the quench to measure the number of soli tons; it must be 
long enough that the soli tons are formed such that they can be 
distinguished from density fluctuations due to noise, but short 
enough that the number of solitons has not appreciably de¬ 
cayed. Fortunately, if the system has periodic boundary con¬ 
ditions, a net winding of the system can emerge post-quench 
as a topologically stable [28] remnant of the domains that can 
also be used to test the theory [10, 23] 

We consider a ID Bose gas in a toroid of circumference 
L. At the freeze-out time the domain size is given by (8), 
and thus the number of domains m N ^ Lj^. The probability 
distribution for the phase of one domain is uniform between 
—n and tt, with variance 

X n n2 2 

( 11 ) 

The variance of the phase difference between two neighbour¬ 
ing domains Afl,- = 6>,+i - 0,- is then cr“(A0,) = 27r^/3, and the 
variance of the accumulated phase around the toroid is 

9 , 27r^ 271^ L 2n^ 

(r\ec) = 2 cT\m = (A - 1)— ^ N— ^ J — . (12) 

The winding is related to the accumulated phase by 'W - 
0c and thus the variance of the winding is given by 


of modes with energy less than a specified cutoff (fcut), and 
the incoherent region (I) which contains the thermalised 
high-energy modes; the incoherent region acts as a thermal 
reservoir to the C-field. The SPGPE has been successfully 
numerically implemented [14, 15] despite its complexities. 
When considering systems that are tightly confined in one or 
more dimensions, the numerical implementation can be made 
more efficient by reducing the dimensionality of the SPGPE 
equations of motion [22]. 

The ID SPGPE is obtained by assuming the low-energy 
fraction of the system affording a SPGPE description is in the 
harmonic oscillator ground states for the two tightly trapped 
dimensions, enabling integrating over these dimensions [22]. 
The transverse dimensions are confined by a parabolic trap 
with harmonic oscillator frequency a>j_ and oscillator length 
flj. = ^IWjmLoZ- The resulting equation of motion is 

{S)hdiff(x, t) - !Px|(i + 7)ii^ — -OiA(Jc, t)dt + hdWy{x, f) 

— iVeix, t)il/dt + ihil/{x, t)dWe{x, f)|. (16) 

where (S) denotes Stratonovich integration. The ID projec¬ 
tion operator Vx implements the energy cutoff ecut in the re¬ 
maining dimension x. The Hamiltonian evolution, with reser¬ 
voir chemical potential as energy reference, is generated by 


cr^i-W) = 



1 L 
6|’ 


(13) 


giving the scaling for the standard deviation of the winding 
distribution 


where 


cr(^ 



(14) 


(15) 


Eollowing a quench, the system may be left to return to a sta¬ 
ble state which may contain some non-zero winding. This 
winding can then be measured experimentally by interferom¬ 
etry [29, 30]. The downside to measuring the winding over 
the defects themselves is that the variance of winding number 
is a higher order moment and thus more trajectories are re¬ 
quired to obtain good statistics. Eurthermore, the power-law 
exponent for the winding is half that for defects, making the 
scaling more susceptible to error. 


£il/{x, f) = ['K(x) -H gi|iA(x, f)l^ - 0 (17) 


where "T/jx) = -h^dll2in+Vcxtix) is the single-particle Hamil¬ 
tonian with external potential yext(x), and gi = 2ha)±as is the 
ID interaction strength with s-wave scattering length a,. The 
reservoir is described by chemical potential (yu), temperature 
(T), and cutoff energy (fcut), and the functions 


80 “ ^ 

—i > -O 

)2 2_J 




dB j=\ 




„2/3ec 


-,hj 


Ve(x, t) - -h J" dx's(x - x')dx' j(x', t), 
ih 

j(x, 0 = [>//dxi/f* - iA*<5x(A] , 

2m 


M r 

27tJ 


e(x) = 

Si(k) = erfcx 

M = 


dke‘'‘''Si(k), 




\ V2 

I6na2 


(87ra^) 


2n-1/2 


gy3(ecui-;/) _ J ’ 


(18a) 

(18b) 

(18c) 

(18d) 

(18e) 

(18f) 


B. SPGPE in one dimension 

The stochastic projected Gross-Pitaevskii equation 
(SPGPE) uses C-field methods to model Bose gases at finite 
temperature [12, 15]. The modes of the system are divided 
into two distinct regions; the coherent region (C) consisting 


where Ads = yj2jih?-IrnksT, fA - IjksT, (A>{z,x,a\ - 
2^0 ^^/('^ + '■1*® L^fch transcendent, and 

erfcx(x) s e'*'erfc(x) (19) 

is the scaled complementary error function. The noise terms 
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are Gaussian, with non-vanishing correlations 

{dWlix, t)dWy(x', t)) - ^ ^ 6(x, x')dt, (20) 
' n 

{dWe(x,t)dWsix',t)) = — —e(x - x')dt, (21) 
n 

where 6{x, x') = Yjnec <Pn{z)<p*„(x') is the 5 function for the C 
region. 

The first two terms on the RHS of (16) give the simple 
growth or number-damping SPGPE defined by (17), (18a), 
and describing particle exchange with the I-region [the PGPE 
is recovered as the special case y s 0]. The final two terms 
on the RHS of (16), are the energy damping terms, describ¬ 
ing energy-exchanging interactions between the C-field and 
the thermal cloud without particle exchange. These terms in¬ 
volve a potential, (18b), and depend upon the divergence of 
the current (18c), with strength determined by the rate func¬ 
tion (18d), scattering kernel (18e), and amplitude (18f). Both 
processes involve an associated noise, (20), (21), and satisfy 
the fluctuation-dissipation theorem. The two reservoir inter¬ 
action processes are shown schematically in Eig. 1 (c). 


C. Two-point correlations of the ID SPGPE 


We first consider the equilibrium system with a constant neg¬ 
ative chemical potential. The coefficient (27) is then time- 
independent, making the time-integrals in (26) trivial. Tak¬ 
ing {(f>* {k, t) (f) (k', t')) from (26) and letting f, t' —> oo, while 
keeping \t — t'\ constant and finite, we then transform back to 
position space to obtain the stationary two-point correlation 
function 


(i/c* (x, Of (x',f'))j = 


ksT „l\x- x' 




-,(1 - ill) 


\t-t’\ 


(28) 


where 
G(x, 0 


4/ 


2^ 


dk e 


-ikx 


1 +k^ 


e^'-^erfc 


/2f-x\ 

\2^) 


+ e^'^^erfc 


2t + X 

lyft 


(29) 


and we have identified the steady state correlation length 





(30) 


Although finding analytic solutions to the full ID SPGPE 
(16) is very difficult, some analytic progress can be made for 
the number-damping SPGPE 

Mf(x, 0 = !^A:|(i + t)(jU - -C)f(x, t)dt -H hdWy(x, o|. (22) 

We assume that we are on the symmetric side of the transition 
ip < 0), allowing us to approximate gilfp a; 0, as the number 
density of the C-field is very low. Transforming to momentum 
space we obtain the equation of motion 

/ ^2^2 \ 

Mf {k, t) - (i H- y) u{t) —-— (j) (k, t) dt + hdU (k, t) (23) 
\ 2m I 


and relaxation time 



(31) 


To verify that this expression has the correct limits, we see 
that G(0, t) - erfc( ^/i), so that using erfc(z) e " I ^Jnz for 
large |z| gives 


(f* (x,Of (x, t'))s 


ksT 

^7r(y - i)\p\\t - t’\lh 


(32) 


confirming that t is the relaxation time. Similarly, G(x, 0) = 
and hence at equal times we recover 


where 


cl>ik,t)^i2n) 


- 1/2 


/- 


dxe '^'f (x, 0 


(24) 


(25) 


is the C-field wave function in momentum space and 

dU (k, t) = {27:)-^!^ J dxe^’^^dWy (x, t) 

is the number damping noise in momentum space. The gen¬ 
eral solution to (23) is 

-I- f (26) 

Jo 


<f * (x, 0 f (x', t))s = <f *(x, t)i//(x', 0) = ^^e ^'1''^, (33) 

confirming ^ as the correlation length. Comparing (30) and 
(31) with (2) and (3) respectively reveals the values of the 
critical exponents for the number-damping SPGPE theory, 
namely the correlation length critical exponent v = 1 /2 and 
the dynamical critical exponent z = 2. We also obtain the 
constants of proportionally for the scaling relations 


^0 


TO 


7j2mpo 

h 

7l^o' 


(34) 

(35) 


where 

1 Ih^k^ \ 

A(k,f) = (/ + y)-(—+ IM 01 J. (27) 


We now turn to the dynamics of the linear quench described by 
(1) for which one may find the equal time correlation function 
on the symmetric side of the transition within the same linear 
approximation used above [9]. The number-damped SPGPE 
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FIG. 2. (colour online) Single trajectories of the C-field during a quench (tq = showing (a) number density and (c) phase for the 

number-damping SPGPE, and (b) number density and (d) phase for the full SPGPE. During initial density growth there are numerous solitons, 
evident as low-density notches with an associated phase jump, that quickly decay leaving a persistent current. The decay of solitons appears 
to be more rapid for the full SPGPE simulations. 


has formal solution (26), for (27) with fi(t) - //of(0- Carrying 
out the integrals in (26), and taking the limit |/ro| ^ oo, we 
find that the equal-time correlation function can be written as 

= (36) 

where 

/(x, y) = ^ f dk e'*^-'erfcx (k^ + , (37) 

•\w J ^ ' 


and fp' = |/t(f)l- In the early part of the quench ^(f) « 

I, and we can use erfc(z) ^ j ■\/7 tz to find the asymptotic 
form for y » 1 


and 


fix,y) 




1 

1 + kP I 


e-My 

y 




2wmt) 


(38) 


(39) 


the adiabatic form expected from (33). In the impulse regime 
(37) is transcendental and must be numerically evaluated. 


however, the dependence upon only |x - x'|/f is a clear signa¬ 
ture of universal scaling in the approach to the critical point. 
Numerically, one finds that in the freeze-out regime /(x,y) 
depends only very weakly on y in the neighbourhood y < 1, 
and the functional form freezes into an impulse regime for 
0 < y < 1. For yu > 0 (1 < y) we must rely on numerical 
simulations of the SPGPE. 


III. SIMULATIONS 
A. SPGPE simulations 

We consider a system with toroidal geometry, where the 
radial extent of the system is much smaller than the circum¬ 
ference of the toroid; the transverse trapping frequency is 
Wj./27r = 200Hz, while the circumference is L = 200aj. 
where aj. is the transverse harmonic oscillator width. The di¬ 
mensionally reduced system is then equivalent to a partially 
degenerate ID Bose gas in a homogeneous trap of length L 
with periodic boundary conditions, embedded in a 3D ther¬ 
mal cloud. We use atomic parameters for *^Rb, giving the 
ID interaction parameter gi = Q.Q139tia>±a±. We chose a 
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FIG. 3. (colour online) The condensate number dynamics for several 
quench times, using (a) the number-damping SPGPE, and (b) the full 
SPGPE. The condensate number was extracted from the numerical 
data using the Penrose-Onsager criterion [see (43)]. 


grid of M - 1024 points, which gives ^ 4 points per corre¬ 
lation length prior to and following the quench (the correla¬ 
tion length takes its smallest value prior to the quench). The 
chemical potential is quenched from ji - -ha>± to fi - fno± 
i.e. /To = The temperature is held constant throughout 
all simulations at T = 0.5 T^, where Tc is the transition tem¬ 
perature for the ideal Bose gas confined to a 3D toroidal trap: 


ticb I N 


(40) 


with a> the geometric mean frequency of the toroid 


2nh 

of = —^0) 
niL^ 


(41) 


f(z) the Riemann zeta function, and N the number of parti¬ 
cles [31]. To determine Tc we have used the Thomas-Fermi 
value of the particle number N - jJ-Llgi corresponding to the 
post-quench parameters, giving a value of ^ 14400. For 
the number-damping SPGPE these parameters give an aver¬ 
age C-field population of A^c ~ 14300 and condensate number 
Nq ^ 12600 in equilibrium post-quench, while the inclusion 
of energy-damping reduces these values to Nc ~ 13100 and 
Nq ^ 10200. 

To obtain a thermalized initial state we use the C-field wave 
function i/r(v) = 0, that is then evolved using the ID SPGPE 
with a high number-damping rate (18a) of y = 1 for 1000 
units of the relaxation time (31) to allow the system to come 
to equilibrium with the thermal cloud in its pre-quench state. 
Eor the quench we set the number-damping rate to a value 
more suitable for dynamics, y = 10^^, and begin the chemi¬ 
cal potential ramp (1). Once the quench is finished (f = tq) 



tquji_ 

FIG. 4. (colour online) The results of our self-similarity algorithm 
with respect to the condensate number for the number-damping 
SPGPE and the full SPGPE. The numerical data for the number¬ 
damping (full) SPGPE is represented by blue (red) points, while the 
green (cyan) line is a least squares fit of a power law to the nu¬ 
merical data. The power law exponents are a = 0.5119± 0.0178 
and a = 0.7145 ± 0.0358 for the number-damping SPGPE and full 
SPGPE respectively. 


we allow the system to evolve for a further 10 units of the 
relaxation time, which we have seen is sufficient for any post- 
quench dynamics to cease. Eor the number-damping SPGPE, 
the energy-damping is not included, so an energy-damping 
rate (18f) of A1 = 0 is used for the entirety of the simulation. 
Eor the full SPGPE, we also include energy-damping terms, 
with rate A4 = y = 10“^ for the entirety of the dynamics. We 
note that y ~ A1 for experimentally relevant parameters [15]. 


B. Results and analysis 

The qualitative behavior observed over the course of the 
quench is illustrated by the density and phase of the C-field 
t), shown for example trajectories of the number-damping 
SPGPE and full SPGPE in Eig. 2; post-quench we see the 
emergence of decaying grey solitons in the density, and even¬ 
tually a stable persistent current is evident in the phase. In or¬ 
der to analyse the SPGPE simulations further we first consider 
the condensate population, a quantity that can be computed 
throughout the quench dynamics in our ID system, allowing 
extraction of the freeze-out time f numerically. 

We find the condensate mode and population using the 
Penrose-Onsager criterion [32], starting from the one-body 
density matrix 


p{x, X, t) = {ijfix, tWix', t)), (42) 


where the angled brackets denote ensemble averaging over 
trajectories; in the classical field regime this is equivalent to 
operator averaging as we may neglect commutators. Solving 
the eigenproblem 



, X, t)4>kix', t) = nk{t)^k{x, t). 


(43) 
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FIG. 5. (colour online) The condensate number over time for sev¬ 
eral quench times using (a) the number-damping SPGPE and, (b) the 
full SPGPE. The time axis has been scaled by the freeze-out time 
as predicted by our self-similarity algorithm for each theory. The 
yellow shaded region indicates the impulse regime where KZM ap¬ 
proximates the system dynamics as frozen. For comparison, the to¬ 
tal C-field population is also shown for one of the quenches (black 
dashed line). 


then gives the system orbitals 0j:(x) and their occupations 
tik, with the largest eigenvalue giving the condensate number 
Nq = sup^ iiji and associated wave function (poix). To construct 
the density matrix we performed 10^ trajectories per quench 
time using both the number-damping ID SPGPE and full ID 
SPGPE. The mean condensate number for several quenches 
is shown in Eig. 3. The curves appear to be of the same 
functional form, but with rescaled time axis. This property 
is known as self-similarity, where rescaling the time axis by a 
particular factor will cause the curves to collapse onto a single 
curve; in this case the rescaling factor is the freeze-out time f. 

We find the freeze-out time from these curves based on their 
self-similarity. A reference condensate number is chosen rel¬ 
atively close to where the condensate number initially grows. 
The time at which this reference is reached should be linearly 


FIG. 6. (colour online) Standard deviation of the final winding for 
various quench times. The blue and red data points are the result 
of simulations of the number-damping SPGPE and the full SPGPE 
respectively. The green line is a a best fit for the number-damping 
SPGPE data in the regime that obeys a power law (14), giving an 
exponent of j8 = 0.1236 ± 0.0098. The cyan line is the equivalent for 
the full SPGPE data, giving an exponent of /I = 0.0966 ± 0.0128. 


proportional to the freeze-out time, and thus obey a power law 

at = a (roTg) , (44) 

where a is a constant of order unity. We then use a power-law 
fit to obtain the exponent a of (6) and thus the value of zv. 
This allows us to determine the actual freeze-out time using 
(5) with the constants (34) and tq (35). 

We use a reference condensate number of No = 100, occur¬ 
ring rapidly after initial condensate growth. Figure 4 shows 
the time at which the reference condensate number is reached 
against the quench time for several different quench times. 
Both theories give results that obey a power law, with expo¬ 
nent (6) differing between the theories. A least-squares-fit to 
the results of simulations using the number-damping SPGPE 
gives or = 0.5119 + 0.0178, consistent with the mean-field pre¬ 
diction a = 0.5. A least-squares-fit for the full SPGPE gives 
a = 0.7145 ± 0.0358, a value that is inconsistent with mean 
field theory. 

In Figure 5 we show the condensate number over time for 
the same quenches as Figure 3, this time with the time axis 
rescaled by the freeze-out time (5) found using the values of 
zv found above and the constants (34) and tq (35). The im¬ 
pulse regime t 6 [-f, i] is highlighted to emphasise that con¬ 
densate growth does not begin until the system has entered 
the adiabatic regime. Rescaling the time axis by f results in 
the curves lying on top of one another for both theories. Note 
that the value for tq used in (44), given by Eq. (35), is essen¬ 
tially exact for the number-damping SPGPE, while for the full 
SPGPE it is only an approximation as the additional damping 
terms change the relaxation time, however the power law ex¬ 
tracted from Fig. 4 is unaltered by this choice. We have also 
included the growth of the C-field particle number Nc(t) for 
one value of the quench time for comparison. 
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FIG. 7. (colour online) The two-point correlation function p{r, t) = (if/* (x, t) if/ (x', t)), where r = x - x', at various times for three different 
quench times simulated using both the number-damping SPGPE and the full SPGPE. The quench times are TgtUj. = (a, d), tqijj± = (b, 

e), and tqld^ = e* (c, f). The top row (a-c) shows the normalised numeric correlation function at ? = -f (blue), f = 0 (red), and t = f (magenta) 
resulting from simulations of the number-damping SPGPE, as well as the normalised analytic correlation function (36) at f = -F (green) and 
f = 0 (cyan). The bottom row (d-f) shows the normalised numeric correlation function at ? = 0 (red) resulting from simulations of the full 
SPGPE, as well as the normalised correlation function (36) at r = 0 (cyan). 


The final winding number is calculated using 'W - 
where Oc is the accumulated phase of the final wave function 
around the toroid. For the purposes of obtaining the statisti¬ 
cal distribution of the final winding number we performed 10"^ 
trajectories per quench time using both the number-damping 
ID SPGPE and the full SPGPE. Eigure 6 shows the stan¬ 
dard deviation of the final winding for a range of quench 
times. Ear from the fast quench limit {tq - 10^) the wind¬ 
ing standard deviation from simulations of both theories obeys 
a power law with respect to the quench time. The exponent 
(15) differs between the two; the number-damping SPGPE 
gives p - 0.1236 + 0.0098, while the full SPGPE gives 
P - 0.0966 + 0.0128. The mean field value/? = 0.125 is within 
error of the number-damping SPGPE value, but not within er¬ 
ror of the full SPGPE. 

The equations (6) and (15) are a pair of simultaneous equa¬ 
tions that relate a and p to the critical exponents v and z, giving 


z 


a _ ip 

1 -a’ 


(45) 


and we can thus determine the values of v and z from our 
numerical simulations of quenches. The number-damping 
SPGPE gives the critical exponents v = 0.5065 ± 0.0586 and 
z = 2.071 + 0.236. These are consistent with the equilibrium 
mean field critical exponents v = 1 /2 and z = 2, a result that 
was also found in Ref. [10]. The full SPGPE give the criti¬ 
cal exponents v = 0.6767 + 0.1745 and z = 3.698 + 0.675, a 
significant departure from the equilibrium mean field theory. 

We have also numerically calculated the two-point cor¬ 
relation function using ensembles of trajectories for several 
quench times; this is shown in Eigure 7. Eor the number¬ 


damping SPGPE, we show the numerical correlation function 
at the times f = {-F, 0, t] with the analytical correlation func¬ 
tion (36) at the times t - {-f, 0}, where t is calculated using 
the constants (34) and to (35) and the numerically obtained 
critical exponents v and z. Eor the full SPGPE, we do not 
have values for the constants and tq, and thus we cannot 
calculate the freeze-out time t. Hence we only show the nu¬ 
merical correlation function at f = 0 and compare this to the 
analytical form for the correlation function (36) at the critical 
point (f = 0). We see that for number-damping at the time 
t - —t and t - 0, the numerical data shows excellent agree¬ 
ment with the analytical expression (36), while at f = F the 
correlation length has clearly grown beyond what is predicted 
by (36); this is unsurprising as (36) was derived under the as¬ 
sumption that the transition has not yet been reached. Eor the 
full SPGPE we see that the correlation length is less than that 
predicted by (36), possibly a consequence of the extra noise. 


IV. DISCUSSION 

Our results indicate that the energy-damping reservoir in¬ 
teraction can have a significant effect on the Bose-Einstein 
condensation transition. While the mean field critical expo¬ 
nents V = 1 /2 and z = 2 were consistent with the number¬ 
damping SPGPE (as can also be shown analytically), the in¬ 
clusion of the energy-damping terms resulted in an increase in 
both V and z to the point where they were no longer consistent 
with mean field theory. Other possible universality classes in¬ 
clude the E-model [33] (v = 2/3, z = 3/2), for which there 
has been experimental evidence [34], or the 3D XY model 
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Exponent 

a = zv/( 1 -1- zv) 

fi = v/2(l ± zv) 

Number-Damping SPGPE 

0.5119 ±0.0178 

0.1236 ±0.0098 

Full SPGPE 

0.7145 ± 0.0358 

0.0966 ± 0.0128 

Mean Field 

0.5 

0.125 

F-Model [33] 

0.5 

0.1667 

3D XY [35] 

0.5677 

0.1452 


TABLE I. The power law exponents for the freeze-out time (5) 
and winding standard deviation (14), for both the number-damping 
SPGPE and the full SPGPE. We have also included the values of 
these exponents as predicted by related universality classes. 


(v = 0.6717, z - 1.9550), which is thought to be the class 
to which the BEC transition belongs [35]. The value of the 
correlation-length critical exponent v from our full SPGPE re¬ 
sults is consistent with both these universality classes, how¬ 
ever the value of the dynamical critical exponent z is not. 

Experimental temperature quenches have been performed 
in toroidally trapped BECs [30] finding power-law exponents 
in agreement with mean field theory. However, a recent ex¬ 
periment extracting the freeze-out correlation length (8) of a 
quenched Bose gas in a box trap [34] produced power-law ex¬ 
ponents more consistent with the E-model [33]. Our results 
suggest that the microscopic reservoir interactions, in partic¬ 
ular the energy damping process, induce non-imiversal mod¬ 
ifications to the critical evolution. The modification is most 
evident in the dynamical critical exponent. Eurther work is 
needed to model specific experiments, and it may pose an ex¬ 
perimental challenge to measure signatures of non-universal 
critical dynamics. In Table I we summarize the power-law ex¬ 
ponents for our simulations and relevant universality classes. 

A further question remains as to the role of the values of 
the two damping rates y and A1; in this work we have as¬ 
sumed they time independent and equal. The precise ratio rel¬ 
evant for experiments, in principle also changing throughout 
the transition, may influence the final prediction for the criti¬ 
cal exponents. A possible extension to this work would be to 
investigate the effects on the critical exponents when varying 
the ratio y/Ad; the true value of this ratio varies depending on 
the reservoir parameters of the experimental system of interest 
(T, 11, fcut) [15]. 

The energy damping terms affecting the condensate popula¬ 
tion may seem surprising given that the energy-damping pro¬ 
cess is number-conserving, however a previous investigation 
has shown that the energy-damping terms provoke a faster, 
more coherent approach to equilibrium [15], as is also con¬ 
sistent with the role of so-called scattering terms in quantum 
kinetic theory [36]. This more efficient equilibration may ex¬ 
plain the lower final winding standard deviation upon inclu¬ 
sion of these terms, as the defects resulting from the phase 
transition are more efficiently damped away. 

The presence of an additional noise source reduces the cor¬ 
relation length at the boundary of the impulse regime (f = —t), 
reducing the domain size and increasing the number of de¬ 
fects. It would hence be informative to investigate the num¬ 
ber of solitons over the course of the quenches, as is done in 


[9], and compare the two theories with experimental data [24]. 
This is a non-trivial task, as distinguishing solitons from den¬ 
sity fluctuations in a noisy system can be difficult, particularly 
soon after the transition when the density is low. Long-time 
evolution to form a stable winding number suggests a coars¬ 
ening dynamics process [37-39]. 

To the best of our knowledge we have given the most com¬ 
plete treatment of the reservoir interactions in the BEC transi¬ 
tion, for the special case of an effectively ID superfluid form¬ 
ing in a ring trap. Despite its microscopic foundation, cer¬ 
tain details of experimental realizations are missing from our 
model. Most notably, our model of the quench ramps the I- 
region chemical potential, but includes no further dynamics 
of the I-region (thus far absent from any SPGPE theory), and 
we do not include a rigorous matching of the system to a spe¬ 
cific set of experimental particle number measurements (as 
has been achieved in systems closer to equilibrium without 
fitting [19], or in a quench by fitting the condensate growth 
rate to experimental data [8]), both of these aspects of mod¬ 
elling experimental quench dynamics remain open problems 
within SPGPE theory. 

As this field is attracting increasing interest, it is of some 
value to connect the SPGPE model to related recent work. We 
note that non-Markovian additive noise was shown to generate 
modified dynamical exponents in a G-L model [40]. Interest¬ 
ingly, quantum quenches may leave important signatures be¬ 
hind in higher order operator moments [41]; such information 
is almost certainly absent from the classical field theory used 
in this work. As the BEC transition is dominated by classical 
fluctuations [42], the SPGPE nevertheless provides a quanti¬ 
tative theory of the transition. 

The comparison of microscopic theories of dissipation with 
generic models also requires some comment. Recent work 
using the holographic duality approach [43] recovered mean- 
field power-laws for the winding number scaling in a model 
of a quenched superconducting transition. While the holo¬ 
graphic method provides a very general approach with certain 
computational advantages, in holographic models the noise 
is introduced by hand via the fluctuation-dissipation theo¬ 
rem (EDT), generating dissipative terms and additive noise 
associated with particle growth. Numerical work typically 
focuses on a regime that is equivalent to overdamped BEC 
dynamics [44], as discussed in Ref. [45]. In contrast, the 
SPGPE theory is derived from a first-principles analysis of 
the Bose gas field theory, includes both additive and multi¬ 
plicative noise (number damping and energy damping), and is 
underdamped; consistency with the EDT is an inherent prop¬ 
erty of the SPGPE. 


V. CONCLUSIONS 

We have simulated chemical potential quenches across the 
Bose-Einstein condensation transition in a ring geometry by 
numerically solving the dimensionally reduced stochastic pro¬ 
jected Gross-Pitaevskii equation with and without the energy 
damping terms. The final winding statistics of the number¬ 
damping SPGPE were found to obey the mean-field power- 
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law exponent predictions of the Kibble-Zurek mechanism. 
The complete SPGPE, including energy damping terms, ex¬ 
hibits a modified power-law for the winding statistics. The 
freeze-out time was also extracted from the condensate num¬ 
ber growth curves, with both theories again resulting in a 
power law with respect to the quench time but with differing 
exponents. While the number-damping SPGPE gave results 
consistent with mean field theory, we were unable to find a 
universality class with critical exponents consistent with the 
complete SPGPE results. The full SPGPE results are strongly 
suggestive of non-universal modifications to critical dynamics 
of the EEC phase transition, as is most clearly demonstrated 


by the modified dynamical exponent z. Our results highlight 
the importance of system-specific reservoir interactions in dy¬ 
namical critical phenomena, as further suggested by recent ex¬ 
perimental studies of the EEC phase transition [24, 30, 34]. 
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